function [Material_State,D_ep]=Elastic_Plastic_Model(Material,Material_State,e)
%Material properties
%--------------------
E=Material.E;
v=Material.v;
f_y=Material.f_y;
%Values from previous iteration
%-----------------
e_i=Material_State.e;
s_i=Material_State.s;
k_i=Material_State.k;
%Constitutive law of elastic material
%----------------------------------------
E_hat=E/((1-2*v)*(1+v));
G=1/2*E/(1+v);
D_e= [E_hat*(1-v) E_hat*v     E_hat*v     0 0 0;...
      E_hat*v     E_hat*(1-v) E_hat*v     0 0 0;... 
      E_hat*v     E_hat*v     E_hat*(1-v) 0 0 0;... 
      0           0           0           G 0 0;...
      0           0           0           0 G 0;...
      0           0           0           0 0 G];
%Strain increment and elastic stress
%-----------------------------------
d_e=e-e_i;
s_e=s_i+D_e*d_e;
%Initial Solution
%------------------
i=0;
s=s_e;
d_k=0;
k=k_i;
while 1
i=i+1;
%Principal Stress
%------------------
[T,S] = eig([s(1) s(4) s(6);...
             s(4) s(2) s(5);...
             s(6) s(5) s(3)]);
S=diag(S);                 
[S,I]=sort(S,'descend'); T=T(:,I);         
%Calculate yield function and its deriviatives
%----------------------------------------------
s_h=f_y+0.1*E*k;   %Hardening Law
f=1/sqrt(2)*sqrt((S(1)-S(2))^2+(S(2)-S(3))^2+(S(3)-S(1))^2)-s_h;
if i==1 && f<0; break; end
TT=[T(1,1)^2 T(2,1)^2 T(3,1)^2 2*T(1,1)*T(2,1) 2*T(2,1)*T(3,1) 2*T(3,1)*T(1,1);...
    T(1,2)^2 T(2,2)^2 T(3,2)^2 2*T(1,2)*T(2,2) 2*T(2,2)*T(3,2) 2*T(3,2)*T(1,2);...
    T(1,3)^2 T(2,3)^2 T(3,3)^2 2*T(1,3)*T(2,3) 2*T(2,3)*T(3,3) 2*T(3,3)*T(1,3)];
    
d_f_d_S= [-(2^(1/2)*(2*S(2) - 4*S(1) + 2*S(3)))/(4*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2));...
          -(2^(1/2)*(2*S(1) - 4*S(2) + 2*S(3)))/(4*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2));...
          -(2^(1/2)*(2*S(1) + 2*S(2) - 4*S(3)))/(4*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2))];...
 
d_d_f_d_d_S=[ 2^(1/2)/((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2) - (2^(1/2)*(2*S(2) - 4*S(1) + 2*S(3))^2)/(8*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(3/2)), - 2^(1/2)/(2*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2)) - (2^(1/2)*(2*S(1) - 4*S(2) + 2*S(3))*(2*S(2) - 4*S(1) + 2*S(3)))/(8*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(3/2)), - 2^(1/2)/(2*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2)) - (2^(1/2)*(2*S(1) + 2*S(2) - 4*S(3))*(2*S(2) - 4*S(1) + 2*S(3)))/(8*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(3/2));...
              - 2^(1/2)/(2*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2)) - (2^(1/2)*(2*S(1) - 4*S(2) + 2*S(3))*(2*S(2) - 4*S(1) + 2*S(3)))/(8*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(3/2)),                          2^(1/2)/((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2) - (2^(1/2)*(2*S(1) - 4*S(2) + 2*S(3))^2)/(8*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(3/2)), - 2^(1/2)/(2*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2)) - (2^(1/2)*(2*S(1) + 2*S(2) - 4*S(3))*(2*S(1) - 4*S(2) + 2*S(3)))/(8*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(3/2));...
              - 2^(1/2)/(2*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2)) - (2^(1/2)*(2*S(1) + 2*S(2) - 4*S(3))*(2*S(2) - 4*S(1) + 2*S(3)))/(8*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(3/2)), - 2^(1/2)/(2*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2)) - (2^(1/2)*(2*S(1) + 2*S(2) - 4*S(3))*(2*S(1) - 4*S(2) + 2*S(3)))/(8*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(3/2)),                          2^(1/2)/((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(1/2) - (2^(1/2)*(2*S(1) + 2*S(2) - 4*S(3))^2)/(8*((S(1) - S(2))^2 + (S(1) - S(3))^2 + (S(2) - S(3))^2)^(3/2))];
 
d_f_d_s=TT'*d_f_d_S;
d_d_f_d_d_s=TT'*d_d_f_d_d_S*TT;
h=0.1*E; %Hardening Law Derivative
d_f_d_k=-h;
h_c=1;
d_h_c_d_S=[0 ; 0 ; 0];
d_h_c_d_s=TT'*d_h_c_d_S;
%Calculate Risidue
%-------------------
r1=D_e\(s-s_e)+d_k*h_c*d_f_d_s;
r2=f;
r=[r1 ; r2];
if norm(r)<10^(-6)*f_y; break; end
%Calculate Jacobian
%--------------------
J11=inv(D_e)+d_k*(d_f_d_s*d_h_c_d_s'+h_c*d_d_f_d_d_s);
J12=h_c*d_f_d_s;
J21=d_f_d_s';
J22=d_f_d_k;
J=[J11 J12 ; J21 J22];
%Calculate Variation of Stress Vector and Hardening Parameter
%--------------------------------------------------------------
dd_sk=-J\r;
dd_s=dd_sk(1:6);
dd_k=dd_sk(7);
s=s+dd_s;
d_k=d_k+dd_k;
k=k_i+d_k;
end
Material_State.e=e;
Material_State.s=s;
Material_State.k=k;
%Calculate Tangential Stiffness Matrix
%---------------------------------------
if i==1 && f<0; D_ep=D_e;
else
H=inv(inv(D_e)+h_c*d_k*d_d_f_d_d_s);
D_ep=H-(h_c*H*(d_f_d_s*d_f_d_s')*H)/(h+h_c*d_f_d_s'*H*d_f_d_s);
end
end
